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Abstract 

We present the latest version of micrOMEGAs, a code that calculates the relic density of the 
lightest supersymmetric particle (LSP) in the minimal supersymmetric standard model (MSSM). 
All tree-level processes for the annihilation of the LSP are included as well as all possible 

'-^ ' coannihilation processes with neutralinos, charginos, sleptons, squarks and gluinos. The cross- 

sections extracted from CalcHEP are calculated exactly using loop-corrected masses and mixings 

^ . as specified in the SUSY Les Houches Accord. Relativistic formulae for the thermal average 

are used and care is taken to handle poles and thresholds by adopting specific integration 
routines. The input parameters can be either the soft SUSY parameters in a general MSSM or 
the parameters of a SUGRA model specified at some high scale (GUT). In the latter case, a link 
with Suspect, SDFTSUSY, Spheno and Isajet allows to calculate the supersymmetric spectrum, 
Higgs masses, as well as mixing matrices. Higher-order corrections to Higgs couplings to quark 
pairs including QCD as well as some SUSY corrections (Am&) are implemented. Routines 
calculating (g — 2) M , b — ► S7 and B s — > u + n~ are also included. In particular the b — ► S7 routine 
includes an improved NLO for the SM and the charged Higgs while the SUSY large tan (3 effects 
beyond leading-order are included. This new version also provides cross-sections for any 2 — * 2 
process as well as partial decay widths for two-body final states in the MSSM allowing for easy 
simulation at colliders. 
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1 Introduction 



We present micrOMEGAsl . 3, a program which calculates the relic density of the lightest 
super symmetric particle (LSP) in the minimal supersymmetric standard model (MSSM). 
The stable LSP, which occurs in supersymmetric models with R parity conservation, 
constitutes a good candidate for cold dark matter. Recent measurements from WMAP[1] 
have in fact constrained the value for the relic density within 10%, 

.094 < Qh 2 < .128 at 2a. 

Forthcoming experiments by the PLANCK satellite [2] will pin-down this important pa- 
rameter to within 2%. One therefore needs to evaluate the relic density with high accuracy. 

The relic density calculation is based on solving the equation characterizing the evo- 
lution of the number density of the LSP. For this, one needs to evaluate the thermally 
averaged cross-section for annihilation of the LSP, as well as, when necessary, coannihila- 
tion with other supersymmetric (SUSY) particles [3, 4, 5]. We use, as in micrOMEGAsl. 1 
[6], the method described in [7] for the relativistic treatment of the thermally averaged 
cross-section, and the generalization of [8] to the case of coannihilation. However we have 
improved our method for solving the density evolution equation, it is now solved numer- 
ically without using the freeze-out approximation. This improvement has not impaired 
the speed of the calculation. 

The other main improvement to micrOMEGAsl . 3 is the use of loop corrected super- 
particle masses and mixing matrices. These masses and mixing matrices, as specified in 
the SUSY Les Houches Accord (SLHA)[9], are then used to compute exactly all annihi- 
lation/coannihilation cross-sections. This can be done whether the input parameters are 
specified at the weak scale or at the GUT scale in the context of SUGRA models or the 
like. In the last case, loop corrections are obtained from one of the public codes which 
calculate the supersymmetric spectrum using renormalization group equations (RGE) 
[10, 11, 12, 13]. These corrections to masses are critical for a precise computation of the 
relic density in two specific regions: the coannihilation region and the region where anni- 
hilation through a Higgs or Z exchange occurs near resonance. Note that these regions of 
the supersymmetric parameter space are among the ones where one predicts sufficiently 
high annihilation rates for the neutralino LSP to meet the WMAP upper bound on the 
relic density. In the first case, the critical parameter is the NLSP-LSP mass difference, in 
the latter the mass difference 2M LS p — m H / z [14]. The Higgs masses are calculated either 
by one of the RGE codes or with FeynHiggsFast[15]. When annihilation occurs near a 
Higgs resonance, higher-order corrections to the width also need to be taken into account. 
As in micrOMEGAsl. 1, QCD corrections to Higgs partial widths are included, furthermore 



we have added the important SUSY corrections, the Amb correction, that are relevant at 
large tan (3. These higher-order corrections also affect directly the Higgs-gg vertices and 
are taken into account in all the relevant annihilation cross-sections. 

Besides the relic density measurement, other direct or indirect precision measurements 
constrain the supersymmetric models. In our package we calculate the supersymmetric 
contribution to (g — 2)^ and to Ap. We also include a new calculation of the supersym- 
metric contribution to B s — > fi + fi~ and an improved calculation of the b — > 57 decay rate. 
The latter includes an improved NLO for the SM and the charged Higgs contribution 
as well as the SUSY large tan/3 effects beyond leading-order, the Amj correction. The 
b — > S7, B s — > fi + fi~ or (g — 2) M routines can be replaced or used as a stand-alone code. 

Within micrOMEGAsl . 3 all (co-)annihilation cross-sections are compiled by CalcHEP [16] 
which is included in the package. CalcHEP is an automatic program to calculate tree-level 
cross-sections for any process in a given model, here the MSSM. We provide a code that 
performs the calculation of cross-sections and decay widths that can be called indepen- 
dently of the relic density calculation. The input parameters are the parameters of the 
SUSY Les Houches Accord. Another new feature is the possibility to call CalcHEP, directly 
from a micrOMEGAsl .3 session, and calculate interactively cross-sections for any process 
in the MSSM or in mSUGRA models. For this, all widths of supersymmetric particles 
are evaluated automatically at tree-level, including the available two-body decay modes. 
The relic density as well as other constraints are also calculated in the CalcHEP session. 

In summary, the new program micrOMEGAsl .3 

• Calculates complete tree-level matrix elements for all subprocesses. 

• Includes all coannihilation channels, in particular channels with neutralinos, charginos, 
sleptons, squarks and gluinos. 

• Calculates the relic density for any LSP, not necessarily the lightest neutralino. 

• *Deals with two sets of input parameters: parameters of the MSSM understood to 
be specified at the EWSB scale or parameters of the SUGRA model specified at the 
GUT scale. Both mSUGRA or non-universal SUGRA models are included. 

• *Provides an interface with the main codes to calculate the supersymmetric spec- 
trum: Suspect [10], Isajet [13], Spheno [12] and SOFTSUSY [11]. 

• *Includes an interface with the SUSY Les Houches Accord [9] for supersymmetric 
model specifications and input parameters. This gives a lot of flexibility as any 
model for which the MSSM spectrum is calculated by an external code can be 
incorporated easily. 



• *Includes loop corrected sparticle masses and mixing matrices. 

• *Includes loop-corrected Higgs masses and widths. QCD corrections to the Higgs 
couplings to fermion pairs are included as well as, via an effective Lagrangian, the 
Anib correction relevant at large tan/5. 

• *Provides exact numerical solution of the Boltzmann equation by Runge-Kutta. 

• Outputs the relative contribution of each channel to 1/Q 

• * Computes cross-sections for any 2-^2 process at the parton level. 

• * Calculates decay widths for all particles at tree-level including all 1 — > 2 decay 



• *Calculates NLO corrections to b — > 57. 

• *Calculates constraints on MSSM: (g — 2) M , Ap, B s — > /i + p~. 

• *Supports both C and Fortran. 

• Performs rapidly the relic density calculation, the limiting factor in the execution 
time of the program is the computation of the supersymmetric spectrum. 

New features in the list above are denoted by a star. In this paper we emphasize mainly 
the new features of our package, full details can be found in the original reference [6]. In 
Section 2, we describe the main changes to our calculation of the relic density. We then 
give the parameters of the supersymmetric model used in our package. A description of the 
package follows in Section 4. Section 5 gives instructions for running the program as well as 
sample sessions. Finally in Section 6 we compare our results with those of DarkSUSY4.0 
[17], the other public package that computes the relic density of supersymmetric dark 
matter. 

2 Calculation of the relic density 

The most complete formulae for the calculation of the abundance Y(T) were presented 
in [7, 8] and we will follow their approach rather closely. The evolution equation for the 
abundance, defined as the number density divided by the entropy density, writes 



modes. 




(2.1) 



where g* is an effective number of degree of freedom [7], M p is the Planck mass and 
Y eq (T) the thermal equilibrium abundance. < av > is the relativistic thermally averaged 
annihilation cross-section of superparticles summed over all channels, 

J29i9j I dsy/sK^^/s/T^aijis) 
< av >= 5 , (2.2) 

i 

where g,i is the number of degree of freedom, a^ the total cross section for annihilation 
of a pair of supersymmetric particles with masses rrii, rrij into some Standard Model 
particles, and Pij{\fs) is the momentum (total energy) of the incoming particles in their 
center-of-mass frame. 

Integrating Eq. 2.1 from T = oo to T = T leads to the present day abundance Y(T ) 
needed in the estimation of the relic density, 

^ = f ^(100(Sl/Mpc)). ^ y ' r °» = 2742 X 10 W 5 ™ < 2 - 3 > 

where s(T ) is the entropy density at present time and h the normalized Hubble constant. 
The present-day energy density is then simply expressed as Plsp — 10.54Qh 2 (GeV/m 3 ). 

Let us rewrite Eq. 2.1 in terms of X = T/M LS p 

dY 

— = A(X)(Y„(Xf - Y(Xfl (2.4) 



A[x) = M^ ^M^IX) ^ < av > p 5) 

First note that one will always have Y(X) Y eq (X) when A(X)Y eq (X) > 1. This is 
the case at X < 1 since the equilibrium abundance Y eq (X) (9(1) [7, 8] and for a typical 
electroweak cross-section , < av >^ lCT^GeV" 2 , and LSP mass, M LS p ~ lOOGeV , one 
has A 10 10 . Choosing a starting point for the solution of the numerical equation at 
small X will rapidly return the solution Y = Y eq . On the other hand when X > 1, Y eq (X) 
decreases exponentially as e~ x . Then neglecting the dependence on X in both A(X) and 
Y eq {X)e x we get 

AY = Y(X) - Y eq {X) = ± (2.6) 

where AY <^ Y eq . In this approximation, AY does not depend on X, whereas Y eq (X) 
decreases exponentially. This can be used to find a starting point X^ for the numerical so- 
lution of the differential equation (2.4). In order to find X^ where AY(XfA = 5 Y eq (XfA 
one can solve 

Y eq (X fl )' = A(X h ) * Y eq (X fl f5(5 + 2) (2.7) 



In the darkOmega routine we use this equation to find X^, Y(Xf 1 ) corresponding to 
8 = 0.1 and solve the differential equation (2.4) by the Runge-Kutta method starting 
from this point [18]. We stop the Runge-Kutta run at point Xf 2 where 

Y eq (X h ) < ±Y(X h ) . (2.8) 



+ / A(X)dX . (2.9) 

J Xt„ 



Then we integrate Eq. 2.4 neglecting the term Y eq (X) 

Y(X )~Y(X f2 )^Jx h 

Note that the temperature T = 2.725/T corresponds to X ss 10 14 . Thus without loss of 
precision we can set X = oo for evaluating Y since A(X) oc 1/X 2 . 

Another routine darkOmegaFO performs the calculation in the freeze-out approxima- 
tion 2 . Here we choose 8 = 1.5 as in Ref.[7] and omit the Runge-Kutta step (Xf = Xf 1 = 
Xf 2 ). The precision of this approximation is about 2% although in some exotic cases the 
approximation works badly. 

As in micrOMEGAsl.l, we include in the thermally averaged cross-section , Eq. 2.2, 
only the contribution of processes for which the Boltzmann suppression factor, B, is 
above some value B f 

K^rrii + m^/T) _ x (m i+mj -2m LSp) 
B = — , ' . fa e m LSP > B e (2.10) 

where m^rrij are the masses of the incoming superparticles. The recommended value is 
B e = 10~ 6 [6]. 

In our program we provide two options to do the integrations, the fast one and the 
accurate one. The fast mode already gives a precision of about 1% which is good enough 
for all practical purposes. The accurate mode should be used only for some checks. 
In the accurate mode the program evaluates all integrals by means of an adaptative 
Simpson program. It automatically detects all singularities of the integrand and checks 
the precision. In the case of the fast mode the accuracy is not checked. We integrate the 
squared matrix elements over the scattering angle by means of a 5 points Gauss formula. 
For integration over s, Eq. 2.2, we use a restricted set of points which depends whether we 
are in the vicinity of a s-channel Higgs/Z/W resonance or not. We increase the number 
of points if the Boltzmann factor corresponding to m po i e is larger than 0.01£> e . 



2.1 Decays of the Higgs scalars 



When the LSP is near a Higgs resonance, it annihilates very efficiently The value of 
the neutralino annihilation cross-section depends on the total width if this width is larger 
2 This function was used in the original version of micrOMEGAs[6]. 



than Rs Tf/Xf, the freeze-out temperature. This is usually the case for large Higgs masses 
of ITeV especialy at large tan (3 due to the enhancement in the bb channel. However the 
width of h(H, A) — > bb receives important QCD corrections. Typically for the heavy 
Higgses (m H > ITeV) the partial width into qq can vary easily by a factor of 2 from the 
tree-level prediction, due mostly to the running of the quark mass at high scale. To take 
these corrections into account we have redefined the vertices hqq, Hqq and Aqq using an 
effective mass that reproduces the radiatively corrected Higgs decays [19]. The effective 
mass at the scale Q writes 

M 2 eff (Q) = M(Q) 2 [l + 5.67a + (35.94 - lMn f )a 2 + (164.14 - n 7 (25.77 - 0.259n/))a 3 

(2.11) 

where a = a s (Q)/n, the scale of the reaction is set to Q — 2m^o, M(Q) and a s (Q) are the 
quark masses and running strong coupling in the MS'-scheme. We use NNLO expressions 
for the strong coupling constants [20] and for the running quark masses [19, 21]. The 
relation between the MS and the pole quark masses are implemented at three-loops 
[20, 21]. These are relevant for the top quark, since we use the pole mass as input 
following the SUSY Les Houches Accord [9]. For b-quark, although m^m^) is the input 
parameter, it is still necessary to compute the pole mass used as an input parameter to 
some of the RGE codes. We set M e ff(Q) = M po i e at scales where the effective mass 
exceeds the value of the pole mass. 

We also take into account the SUSY-QCD corrections [22] to h,H,A — > bb vertices 
that are important at large tan (3. Here we use the effective Lagrangian 



C e ff = V 4:7Ca QED 



m b 



1 + Am b 2M W sin 9 W 



TTT -cosa / Am b tano; 
-Hbb 1 + 



+iMtmP I - J^U + hbb J- f 1 - 



cos/3 \ tan/5 J 

1 ' (2-12) 



tan (3 2 J cos (3 \ tan a tan /^ 

where is the effective b-quark mass described above, oiqed the electromagnetic cou- 
pling, tan (3 is the ratio of the vev's of the Higgs doublets and a is the Higgs mixing angle. 
Amj is a correction factor arising from loop contribution of SUSY particles. This factor 
is particularly important at large tan/3 and also contributes to b — > 57 (all details are 
given in Appendix B). 

In the large tan (3 case, when neutralino annihilation via s-channel Higgs exchange 
dominates, the inclusion of SUSY-QCD corrections can shift by about 15% the value for 
the relic density. There is an option to switch off this correction (see Section4.1). 

The total width of the Higgs includes only the two-body final states that occur at 
tree-level. In the case of the light Higgs, this underestimates the width since the partial 



width to off-shell W or gg final states can reach 10%. However an accurate value for this 
very narrow width has in general not a strong impact on the relic density. On the other 
hand a precise value for the heavy Higgs width is necessary. 

2.2 Neutralino "width" 

We assume that the LSP is stable because of R-parity conservation, however it is necessary 
to introduce a width for this stable particle in order to avoid infinities in some processes. 
For example, in the coannihilation process like eiXi ~^ e ^ Yia - t-channel exchange of \i 
an infinity is caused by the pole in the propagator, this is due to the fact that one can 
have a real decay e L — > ex±- We assign a value of sWidth ■ M LSP to the width of all 
super symmetric particles . The default value for the variable sWidth is 0.01. 

2.3 Loop corrections to the MSSM spectrum. 

In the mSUGRA model, but also in the more general MSSM, annihilation of the LSPs near 
a Higgs or Z resonance and/or coannihilation processes are often the dominant reactions 
in models where Vth 2 as 0.1 [23]. For an accurate calculation of the relic density it is then 
very important to have the exact relations between particle masses. In particular, the 
direct annihilation of a pair of neutralinos (x?) depends sensitively on the mass difference 
with the Higgs or Z, 2m^o — M H / Z , when the annihilation occurs near the resonance. 
Furthermore coannihilation processes depend strongly on the NLSP-LSP mass difference. 

In this new version of micrOMEGAsl . 3, we provide an option to calculate loop cor- 
rections to all sparticle masses 3 . Within the MSSM defined at the EWSB scale, loop 
corrections are implemented by a call to Suspect [10], within the SUGRA or other model 
defined at the GUT scale, the loop corrections are done by any of the four public codes 
(Suspect ,S0FTSUSY , Spheno ,Isajet ) for calculating the supersymmetric spectrum 
based on renormalization group equations. Because it is a mass difference rather than 
the absolute mass that has a large impact on the prediction of the relic density, even 
radiative corrections at the percent level, such as is often the case for neutralinos, need to 
be taken into account. Indeed large shifts in the prediction of the relic density between 
tree-level and loop-corrected masses can be found. Typically the prediction for the relic 
density can change by 20%, but in some scenarios corrections can reach 100% or even 
more. We not only use the loop-corrected sparticle masses but also the corresponding 
mixing matrix elements. In this way we take into account some of the loop corrections in 
the evaluation of the matrix elements for different processes. This however means, since 
3 Pole masses in the calculation of the relic density were first used in Ref. [24] 



Table 1: Standard Model parameters 



name 


default 


definition 


AlfEMZ 


0.00781653 


electromagnetic coupling a em (M z ) 


AlfSMZ 


0.1172 


strong coupling, ocf (Mz) for rif — 5 


SW 


0.481 


Weinberg angle, sm9w 


MZ 


91.1884 


Z mass 


Ml 


1.777 


tau-lepton pole mass 


Mtp 


175.0 


t-quark pole mass 


MbMb 


4.23 


MS scale independent b-quark mass Mb (Mb) 



it is only a partial implementation of loop corrections, that theoretical inconsistencies in 
the model could occur, in particular problems with unitarity violation in some processes. 
This would mainly show up in processes with production of gauge particles, however at 
much higher energies that are typically involved in the LSP annihilation processes. 

3 The MSSM parameters. 

In our package, we compute various matrix elements and cross-sections for 2 — > 2 pro- 
cesses within the framework of the MSSM. The model file corresponding to the specific 
implementation of the MSSM was derived with LanHEP[25], a program that generates 
the complete set of particles and vertices once given a Lagrangian [26, 27]. Names are 
attached to the parameters of the MSSM, including those of the SM, and their values can 
be set with an instruction. For example, the command assignVal ( "Mtp" , 180 . ) assigns 
the value m t = 180GeV to the pole mass of the t-quark. 

The list of parameters of the Standard Model and their default values is presented in 
Table 1. All quarks and leptons of the first two generations are assumed massless. The 
default values for the electromagnetic coupling and the Weinberg angle correspond to the 
values in the MS scheme at the M z scale. 

The parameters of the MSSM are described in Table 2. We follow the conventions 
of the SUSY Les Houches Accord [9]. The masses of the third generation fermions are 
ordered, for example corresponds to the lightest top-squark. In this list, the number 
of parameters exceeds the number of MSSM independent parameters. They correspond 
to physical parameters, masses and mixings. This extended set of parameters is however 
necessary when one wants to use effective masses and vertices that include loop corrections. 
Our computation of matrix elements for cross-sections is based on this set of parameters. 
Note that the trilinear muon coupling, A^, is added to the parameter list even though 
it does not contribute to matrix elements or to the spectrum since the muon is assumed 



to be massless. This parameter is however important for evaluating the muon anomalous 
magnetic moment. 



Table 2: MSSM parameters of the SUSY Les Houches Accord 



name 


comment 


name 


comment 


tb 


tan/3 


MSnl 


r-sneutrino mass 


alpha 


Higgs a angle 


MSe^ 


masses of left / right selectrons 


mu 


Higgs \x parameter 


MSm| 


left /right smuon masses 


Mh 


Mass of light Higgs 


MSli 


i= 1,2 masses of light/heavy f 


MH3 


Mass of CP-odd Higgs 


MSu^ 
MSsf 


masses of left / right u-squarks 


MHH 


Mass of Heavy Higgs 


masses of left / right s-squarks 


MHc 


Mass of charged Higgs 


MSti 


i=l,2 masses of light /heavy t-squarks 


Al 


f trilinear coupling 


MSd| 
MScf 


masses of left / right d-squarks 


Am 


jl trilinear coupling 


masses of left / right c-squarks 


Ab 


b trilinear coupling 


MSbi 


i= 1,2 masses of light/heavy b-squarks 


At 


t trilinear coupling 


Zn i:j 


i,j=l,..,4; neutralino mixing matrix 


MNEi 


i=l,2,3,4; neutralino masses 


Zuij 


i=l,2;j=l,2; chargino U mixing matrix 


MCi 


i= 1,2 chargino masses 


Zv i:j 


i=l,2;j=l,2; chargino V mixing matrix 


MSG 


mass of gluino 


Zly 


i=l,2;j=l,2; f mixing matrix 


MSne 


e-sneutrino mass 


Ztjj 


i=l,2;j=l,2; t mixing matrix 


MSnm 


/x-sneutrino mass 


Zbij 


i=l,2;j=l,2; b mixing matrix 



The values of the SLHA parameters can either be set by an external program, here a 
call to one of the RGE codes that calculate the supersymmetric spectrum, or by specifying 
the MSSM parameters at the weak scale. In either case one needs to specify a set of 
independent parameters as described below. 

3.1 Input parameters at the GUT scale 

Within the context of the SUGRA scenario for supersymmetry breaking the MSSM pa- 
rameters can be evaluated at the weak scale starting from a set of scalar masses, gaugino 
masses, trilinear couplings defined at the GUT scale. The GUT scale input parameters 
are listed in Table 3. Only one parameter, tan/5, is defined at Mz- We implicitly as- 
sume that the first two generations are identical. The parameters for the mass of the 
Higgs doublet can be entered with a negative sign, in this case they will be understood 
as M 2 Hu = -\M Hu \ 2 

We treat the mSUGRA model as a special case of the general SUGRA. Since simplify- 
ing relations are imposed on masses and couplings, in the mSUGRA model one has to spec- 
ify only a small number of input parameters at the GUT scale: M , My 2 , A , tan (3, sgn{mu) 
These correspond to 



mO = Mli = Mri = Mqi = Mui = Mdi = MHu = MHd 

- common scalar mass at GUT scale; 
mhf = MG1 = MG2 = MG3 - common gaugino mass at GUT scale; 
aO=At = Ab = Al - trilinear soft breaking parameter at GUT scale; 
tb- tan/? or the ratio of vacuum expectation values at MZ; 
sgn- +/-1, sign of /i, the Higgsino mass term. 

Four different routines read the parameters of Table 2 and pass them to the corre- 
sponding packages that solves the RGE equations and calculate the MSSM masses and 
mixing matrices. The routines suspectSUGRA [10],sof tsusySUGRA [11], sphenoSUGRA [12], 
isajetSUGRA [13] are described in section4.1. Note that some of the standard parameters 
of Table 1 also play a role in the low energy boundary conditions implemented in the 
RGE codes. They are passed to RGE routines implicitly. We assume that the second 
generation is identical to the first one and only parameters of the first generation are used. 



Table 3: Independent GUT-scale parameters. 



name 


comment 


name 


comment 


tb 


tan (3 (at M z ) 


Mil 


Left-handed slepton mass for i st /2 nd gen. 


At 


t trilinear coupling 


M13 


Left-handed slepton mass for 3 rd gen. 


Ab 


b trilinear coupling 


Mrl 


Right-handed slepton mass for I st /2 nd gen. 


Al 


f trilinear coupling 


Mr3 


Left-handed slepton mass for 3 rd gen. 


MG1 


U(l) Gaugino mass 


Mql 


Left-handed squark mass for \ st /2 nd gen. 


MG2 


SU(2) Gaugino mass 


Mq3 


Left-handed squark mass for 3 rd gen. 


MG3 


SU(3) Gaugino mass 


Mul 


Right-handed u-squark mass for \ st /2 nd gen. 


sgn 


sign of fi at the EWSB scale 


Mu3 


Right-handed u-squark mass for 3 rd gen. 


MHu 


Mass of first Higgs doublet 


Mdl 


Right-handed d-squark mass for \ st /2 nd gen. 


MHd 


Mass of second Higgs doublet 


Md3 


Right-handed d-squark mass for 3 rd gen. 



3.2 Input parameters at the weak scale 

The parameters of the SUSY Les Houches Accord can also be calculated starting from 
the set of independent MSSM parameters at the EWSB scale 4 listed in Table 3 [26]. This 
can be done either at tree-level or with loop corrections (see Section 4.1). The names of 
the independent parameters of the MSSM are identical to the GUT scale parameters safe 
for MHu, MHd which are conveniently replaced by fi and M^(MH3). Furthermore at the 
EWSB scale one must define the sfermion masses for all three generations. Here MH3 and 
MG3 are the pole masses of the CP-odd Higgs and of the gluino. All other parameters 
4 This set of parameters was used in the previous version of micrOMEGAs[6]. 



are treated as running ones. When evaluating loop corrections to pole masses starting 
from the independent set of parameters, it is assumed that the parameters are specified 
in the DR scheme at the EWSB scale, Q = y/m^ ■ m t ~ 2 . 



Table 4: Set of independent MSSM parameters at the weak scale. 



name 


comment 


name 


comment 


tb 


tan (3 


MG3 


SU(3) Gaugino mass (gluino 


mass) 


mu 


Higgs \i parameter 


Mli 


Left-handed slepton mass for 


i th 


generation 


At 


t trilinear coupling 


Mri 


Right-handed selectron mass 


for 


jth g enera tion 


Ab 


b trilinear coupling 


Mqi 


Left-handed squark mass for 


jth 


generation 


Al 


f trilinear coupling 


Mui 


Right-handed u-squark mass 


for 


jth generation 


Am 


jl trilinear coupling 


Mdi 


Right-handed d-squark mass 


for 


jth generation 


MG1 


U(l) Gaugino mass 


MH3 


Mass of Pseudoscalar Higgs 






MG2 


SU(2) Gaugino mass 











Two options are available to specify the weak scale MSSM parameters, either from a 
file using the function ewsblnitFile or directly as argument of the function ewsbMSSM. 
Either option will evaluate the supersymmetric spectrum at tree-level or to one-loop 
according to the value of the parameter LCOn, see section 4.1. 

After evaluation of the spectrum in the context of the SUGRA or MSSM models, the 
function calcDep chooses the lightest supersymmetric particles and calculates the running 
masses of quarks at the LSP scale as well as various widths. 

4 Functions of micrOMEGAs 

The routines presented below belong to the micromegas.a library. They are available 
both in the C and Fortran versions. If for some reason a Fortran call differs from the 
C one, we present the Fortran version in brackets "[ ]". The types of the functions and 
their arguments are specified in Appendix A. Examples of implementation are presented 
in Section 5.4. Note that after assignments of the MSSM parameters the user has to call 
the initialization procedure calcDep (Sec. 4.1). Other routines of the package can only 
be used after making this call. 

4.1 Variable assignment and spectrum calculation 

• as s i gnVal (name , val ) 

changes values of the parameters, name is one of the names presented in Tables 1,2, val 
is the value to be assigned. The function returns when it successfully recognizes the 



parameter name and 1 otherwise, 
•assign ValW (name , val) 

the same routine as assignVal, instead of returning an error code it writes a warning on 
the screen. 

•suspectSUGRA ( tb , MG1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , Mil , M13 , Mr 1 , Mr3 , Mql , Mq3 , 
Mul,Mu3,Mdl,Md3) 

calculates the values of the MSSM parameters in the SUGRA scenario using the Suspect 
package. Returns when the spectrum is computed succesfully, 1 in case of non-fatal 
problems (see the manual for the meaning of non-fatal errors [10]), and (—1) if no solu- 
tion to RGE can be found for a given set of boundary conditions. This routine assigns 
values for the parameters in Table 2. The result depends on the input values of the 
SM parameters, in particular on the quark masses, mf°' 6 , m^rrib) (Mtp, MbMb) and on 
the strong coupling constant a s (M z ) (Alf SMZ). These parameters play a role in the low 
energy boundary conditions and are passed implicitly. 

•sof t SusySUGRA (tb , MG1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , Ml 1 , M13 , Mr 1 , Mr3 , Mql , Mq3 , 

Mul,Mu3,Mdl,Md3) 
same as above for SOFTSUSY. 

•sphenoSUGRA ( tb , MG1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , Ml 1 , M13 , Mr 1 , Mr3 , Mql , Mq3 , 
Mul,Mu3,Mdl,Md3) 

same as above for Spheno. 

• isaj et SUGRA (tb , MG1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , Ml 1 , M13 , Mr 1 , Mr3 , Mql , Mq3 , 

Mul,Mu3,Mdl,Md3) 

same as above for Isajet. This function depends only on mf°' e , other SM parameters, and 
in particular mb(mb) and a s , are fixed internally. Isajet does not calculate the trilinear 
muon coupling, we use the approximate relation for mSUGRA models, = A — 0.7M 1 / 2 . 

Note that only the Suspect code is included in our package. Other codes should be 
installed independently by the user and linked to micrOMEGAs as explained in Section 5.1. 

• ewsbMSSM (tb , MG1 , MG2 , MG3 , Am , Al , At , Ab , MH3 , mu , Ml 1 , M12 , M13 , Mr 1 , Mr2 , Mr3 , 

Mql , Mq2 , Mq3 , Mul , Mu2 , Mu3 , Mdl , Md2 , Md3 , LCOn) 

calculates the supersymmetric spectrum at tree-level or one-loop from the set of indepen- 
dent MSSM parameters at the EWSB scale as specified by the parameter LCOn. The Higgs 
sector parameters, masses and mixing angle a, are calculated with FeynHiggsFast [15]. 

LC0n=0 - tree level formulae for super particles masses; 

LC0n=l - Suspect is used to evaluate loop corrections to masses of super particles. 
•ewsblnitFile (filename , LCOn) 



reads the input file filename which specifies the set of independent MSSM parameters 
at the EWSB scale and calculates the supersymmetric spectrum at tree-level or one-loop 
as set by the parameter LCOn (same as above). 

The function returns: 

- when the input has been read correctly; 
-1 - if the file does not exist or can not be opened for reading; 
-2 - if some parameter from Table 4 is missing as displayed on the sceen; 
-3 - if the spectrum cannot be calculated; 
n - when the line number n has been written in the wrong format. 
For example, the correct format of a line is 
MG3 1500. 

•readLesH(f ilename ,LE) 

reads the input file in the SUSY Les Houches Accord format [9]. If LE=1 the SM parameters 
of Table 1 as well as tan (5 are also read from a SLHA output file. 

•calcDep(dMbOn) 

initializes internal parameters for subsequent calculations. In particular, the running 
masses of quarks, the strong coupling constant as well as the widths of gauge bosons, 
Higgses and superparticles. Running parameters are evaluated at the LSP scale. This 
routine also sorts the superparticles and selects the LSP. The parameter dMb0n= switches 
off SUSY-QCD corrections, Am b , see Section 2.1. 

4.2 Display of parameters. 

•f indVal (name ,&val) [f indVal (name , val) ] 

assigns to the variable val the value of the parameter name. It returns zero if such variable 
indeed exists and 1 otherwise. This function can be applied to any of the parameters in 
Table 1,2 as well as to particle masses and widths specified in Tables 5,6,7. 
•f indValW (name) 

returns the value corresponding to the variable name. If name is not defined f indVal W 
writes a warning on the screen. 
•printVar(f ile.N) [ print Var(N) ] 

prints the first N records of the full list of model parameters. The first 7 parameters 
correspond to Table 1, the following 75 parameters correspond to the list in Table 2. To 
see the parameters on the screen, substitute f ile=stdout. In the Fortran version, only 
display on the sceen is possible. 

•printMasses (file , sort) [ printMasses (sort) ] 



prints into the file the masses of the supersymmetric particles as well as all Higgs masses 
and widths. The Fortran version writes down on the screen. If sort^ 0, the masses are 
sorted in increasing order. 
•IspO [ lsp(name) ] 

returns the name of the LSP. The relic density can be calculated with any particle being 
the LSP even though only the neutralino and the sneutrino can be dark matter candidates. 
If the user wants to impose a specific LSP, the nature of the LSP must be checked after 
calling calcDep. 

•lspmass_() [ IspmassO ] 

returns the mass of the lightest supersymmetric particle in GeV. 

4.3 Calculation of relic density. 

•darkOmega(&Xf, fast, Beps) [ darkOmega(Xf ,Fast,Beps) ] 

This is the basic function of the package which returns the relic density Vth 2 (Eq. 2.3). 
The procedure for solving the evolution equation using Runge-Kutta was described in 
Section 2. The value of the freeze-out parameter Xf is returned by the function and 
equals (X/ x +X/J/2, (see the definition in Eq. 2.7, 2.8). The parameter Beps defines the 
criteria for including a given channel into the sum for the calculation of the thermally 
averaged cross-section, Eq. 2.10; 1CT 6 is the recommended value. 

If fast=0, we use an integration routine that increases the number of points until 
an accuracy of 10~ 3 is reached. If fast=l the accuracy is not checked, but a set of 
points is chosen according to the behaviour of the integrand: poles, thresholds, Boltzman 
suppression at large energy. The accuracy of this mode is about 1%. Finally, fast=2 
corresponds to the calculation of relic density using the widely-used approximation [5] 
based on the expansion in terms of velocity 

p ■ a(p) = A + B ■ p 2 . 

The recommended mode is f ast=l. 

If some problem is encountered, darkOmega returns (—1). 

•dark0megaF0(&Xf, fast, Beps) [ dark0megaF0(Xf , fast , Beps) ] 

calculates the relic density as the function darkOmega described above, but using the 
freeze-out approximation. 

•printChannels (Xf , cut, Beps, prcnt,f) [ printChannels (Xf , cut , Beps ,prcnt) ] 

prints the relative contribution to Q^ 1 for all subprocesses for which this contribution 
exceeds the value chosen for cut. If prcnt=l the contribution is given in percent, otherwise 



the absolute value is displayed. It is assumed that the Xf parameter was first evaluated 
by darkOmega. In the C version, the output is directed to the file f , the Fortran version 
writes on the screen. Actually this routine evaluates the partial contributions to the 
integral of Eq. 2.9 without the l/Yf term and returns the corresponding value for VLh 2 . 

4.4 Routines for constraints. 
•deltarho.O [ delrhoO] 

calculates, by a call to a Suspect routine, the Ap parameter which describes the MSSM 
corrections to electroweak observables. It contains stop/sbottom contributions, as well 
as the two-loop QCD corrections due to gluon exchange and the correction due to gluino 
exchange in the heavy gluino limit [28]. Precise measurements of SM electroweak observ- 
ables allow to set the limit Ap < 2 ■ 1CT 3 . 

•bsgnlo_() [ bsgnloO ] 

returns the value of the branching ratio for b — > sj. For b — > 57 we have improved on the 
results of [29] by including some very recent new contributions beyond the leading order 
that are especially important for high tan/3. Full details can be found in Appendix B. 

•bsmumu_() [ bsmumuO ] 

returns the MSSM contribution to B s — > fi + fi~. Our calculation is based on [30] and 
agrees with [33]. It includes the loop contributions due to chargino, sneutrino, stop and 
Higgs exchange. The Amj effect relevant for high tan (3 is taken into account. The current 
bound from CDF experiment at Fermilab is B.R.(S S — > fi + fi" < 9 x 10~ 6 ) [31] and the 
expected bound from Runlla should reach B.R.(i? s — > fi + fi~ < 2 x 1CT 7 ) [32]. 

•gmuon_() [ gmuonO ] 

returns the value of the supersymmetric contribution to the anomalous magnetic moment 
of the muon [34]. The result depends only on the parameters of the chargino/neutralino 
sector as well as on the smuon parameters, in particular the trilinear coupling (Am). 
Our formulas agree with [35]. The latest experimental data on the (g — 2) M measurement 
using /i _ [36], brings the average to a^ xp ' = 11659208 ±6 x 10~ 10 . The quantity a M includes 
both electroweak and hadronic contributions and is still subject to large theoretical errors, 
the allowed range for 5a^ = a^ xp - — a^ hco - then has also large errors. We estimate the 3a 
range to be 5.1 < Sa^ x 10 10 < 64.1 [37]. 

•masslimits_() [ masslimitsO ] 

returns a positive value and prints a WARNING when the choice of parameters conflicts 
with a direct accelerator limits on sparticle masses. The constraint on the light Higgs 
mass is not implemented and must be added by the user. 



Among the routines that calculate constraints, only masslimits issues a warning if the 
chosen model gives a value outside the experimentally allowed range. All other constraints 
must be checked by the user. 

4.5 QCD auxiliary routines. 

• alphaQCD(Q) 

calculates the running a s at the scale Q in the MS scheme. The calculation is done using 
the NNLO formula in [20]. Thresholds for b-quark and t-quark are included in n/ at the 
scales mfe(mfe) and m t (m t ) respectively. Implicit input parameters are AlfSMZ, Mtp, and 
MbMb defined in Table 1. 

• MbRun(Q) , MtRun(Q) 

calculates top and bottom running masses evaluated at NNLO. 

• MbEff(Q), MtEff(Q) 

calculates effective t- and b-quark masses as in Eq. 2.11. 
•deltaMbO 

calculates the SUSY corrections to Am b (Appendix B). 

4.6 Partial widths and cross sections 

• decay2(pName,k, outl, out2) 

calculates the decay widths (in GeV) for any 1 — > 2 processes. The input parameters are 
pName, the name of the decaying particle and k, the channel number, outl and out2 are 
the names of outgoing particles for channel k. If k exceeds the total number of channels, 
then outl and out2 are filled as empty strings. 

• newProcess (procName , libName) [ newProcess(procName, libName, address) ] 

prepares and compiles the codes for any 2 — > 2 reaction in the MSSM. The result of the 
compilation is stored in the library 

source/ '2-2/ 'libName . os. 
If this library already exists, it is not recompiled and the correspondence between the con- 
tents of the library and the procName parameter is not checked. libName is also attached 
to the names of routines in the UbName.so library. Therefore libName should not contain 
symbols such as +,—,*,/, which are not legal as identifiers. Library names should not 
start with omglib, these are reserved for the libraries used to evaluate Qti 2 . 

The process should be specified in CalcHEP notations, for example 
"e,E->~l+,~l-" 



Table 5: Higgs particles. 



Name 


symbol 


mass 


width 


Name 


symbol 


mass 


width 


Light Higgs 


h 


Mh 


wh 


CP-odd Higgs 


H3 


MH3 


wH3 


Heavy higgs 


H 


MHH 


wHh 


Charged Higgs 


H+,H- 


MHc 


wHc 



without any blank space. One can find all symbols for MSSM particles in Tables 5,6,7. 
Multi-process generation is also possible using the command 

"e,E->2*x" 
where x means arbitrary final states. 

The newProcess routine returns the address of the static structure with contains, for 
further use, the code for the processes. If the process can not be compiled, then a NULL 
address is returned (address[l]=0 in Fortran). newProcess can also return the address 
of a library that was already generated, for example, newProcessC" , "omglib_ol_ol") 
returns the address of the library for neutralino annihilation. 

• inf or22 (address , nsub , nl , n2 , n3 , n4 , &ml , &m2 , &m3 , &m4) 
[ infor22(address,nsub,nl,n2,n3,n4,ml,m2,m3,m4) ] 

allows to check the contents of the library produced by newProcess. Here address is 
the returned value of newProcess call and nsub the subprocess number. The parameters 
returned correspond to the names of particles for a given subprocess (nl, n2, n3, n\) as 
well as their masses (ml, m2, m3, m\). The function returns 2 if the nsub parameters 
exceed the limits and otherwise. 

• cs22 (address , nsub, P, cl, c2 , &err) 

evaluates the cross section for a given 2^2 process with center of mass momentum 
P(GeV). The differential cross section is integrated from cl < cos(6) < cl and 9 is the 
angle between pi,p3 in the center-of-mass frame. If nsub exceeds the maximum value for 
the number of subprocesses then err contains a non zero error code. 

5 Work with the micrOMEGAs package. 

5.1 Installation and link with RGE packages. 

micrOMEGAs can be obtained at 

http : //wwwlapp . in2p3 . f r/lapth/micromegas 

The name of the file downlaoded should be micromegas_l .3.0. tar . gz. After unpacking 
the file, the root directory of the package, micromegas_l . 3 . 0, will be created. This 
directory contains the micro_make file, some sample main programs, a directory for the 



Table 6: Names, masses and widths of supersymmetric particles. 



Name 


symbols 


mass 


width 


Name 


symbols 


mass 


width 


chargino 1 


~1+ 


~1- 


MCI 


wCl 


mu-sneutrino 


nm 




MSnm 


wSnm 


chargino 2 


~2+ 


~2- 


MC2 


wC2 


tau-sneutrino 


~nl 




MSnl 


wSnl 


neutralino 1 


~ol 




MNE1 


wNEl 


u-squark L 


~uL 


~UL 


MSuL 


wSuL 


neutralino 2 


~o2 




MNE2 


wNE2 


u-squark R 


~uR 


~UR 


MSuR 


wSuR 


neutralino 3 


~o3 




MNE3 


wNE3 


c-squark L 


~cL 


~CL 


MScL 


wScL 


j 1 • a 

neutralino 4 


~o4 




TV T~\ TT — 1 A 

MNE4 


wNE4 


c-squark R 


~cR 


~CR 


MScR 


wScR 


gluino 


~g 




MSG 


wSG 


t-squark 1 


~tl 


~T1 


MStl 


wStl 


selectron L 


~eL 


~EL 


MSeL 


wSeL 


t-squark 2 


~t2 


~T2 


MSt2 


wSt2 


selectron R 


~eR 


~ER 


MSeR 


wSeR 


d-squark L 


~dL 


~DL 


MSdL 


wSdL 


smuon L 


~mL 


~ML 


MSmL 


wSmL 


d-squark R 


~dR 


~DR 


MSdR 


wSdR 


smuon R 


~mR 


~MR 


MSmR 


wSmR 


s-squark L 


~sL 


~SL 


MSsL 


wSsL 


stau 1 


~11 


~L1 


MS11 


wSll 


s-squark R 


~sR 


~SR 


MSsR 


wSsR 


stau 2 


~12 


~L2 


MS12 


wS12 


b-squark 1 


~bl 


~B1 


MSbl 


wSbl 


e-sneutrino 


~ne 




MSne 


wSne 


b-squark 2 


~b2 


~B2 


MSb2 


wSb2 



Table 7: Designations for the Standard Model particles 



Name 


symbols 


Mass 


Width 


Name 


symbols 


Mass 


Width 


photon 


A 










tau-neutrino 


nl 


Nl 








Z boson 


Z 




MZ 


wZ 


tau-lepton 


1 


L 


Ml 





W boson 


W+ 


W- 


MW 


wW 


s-quark 


s 


S 








gluon 


G 










c- quark 


c 


C 








electron 


e 


E 








u-quark 


u 


u 








muon 


m 


M 








d-quark 


d 


D 








e-neutrino 


ne 


Ne 








t-quark 


t 


T 


Mt 


wt 


mu-neutrino 


nm 


Nm 








b-quark 


b 


B 


Mb 






source code, a directory for CalcHEP interactive sessions and a directory containing data 
files. To compile, type either 

. /micro_make 

This command is a Unix script, which detects the operating system and its version, 
sets the corresponding compiler options, and compiles the code. Being launched without 
arguments, micro_make compiles only auxiliary libraries needed for relic density evalu- 
ation. Otherwise, the first argument is treated as a C or Fortran main program which 
should be compiled and linked with these libraries. The executable file created has the 
same name as the main program without the .c/.f extension. 

It is interesting to investigate the relic density in the framework of some scenario of 
supersymmetry breaking. We rely on the public codes that evaluate the supersymmetric 
spectrum in the context of models defined at the GUT scale such as the mSUGRA model. 



One of these packages, Suspect [10], is included into the micrOMEGAs package. We also 
support an interface with SOFTSUSY [11], Spheno [12] and Isajet [13]. 

To use Isajet, the corresponding library should be attached to the code. It can be 
done via the variable EXTLIB to be defined in the micro_make file. For example, to use 
Isajet located in the ~/isajet769 directory the definition should be 

EXTLIB="$H0ME/isajet769/libisajet.a" 
If mathlib from CERNLIB is not included in libisajet.a it should be specified in EXTLIB, 
for example 

EXTLIB="$H0ME/isajet769/libisajet .a -L/cern/pro/lib -lmathlib" 

The interface with SOFTSUSY and Spheno is realized in the framework of the SUSY 
Les Houches accord[9] by direct execution of the corresponding programs. In both cases, 
the user has to define in the micro_make file, the variables SOFTSUSY or SPHENO which 
identifies the directory where the corresponding executable file is located. For example, 
S0FTSUSY=$H0ME/sof tsusy.l . 8 

or 

SPHEN0=$H0ME/SPheno2 .2.0 

To install the package, one needs initially about 20MB of disk space. As the pro- 
gram generates libraries for annihilation processes only at the time they are required, the 
total disk space necessary can double after running the program for different models as 
described in the next section. 

5.2 Dynamic generation of matrix elements and their loading. 

In order to take into account all possible processes of annihilation of superparticles into 
SM particles, we need matrix elements for about 2800 different subprocesses. However, 
for a given set of parameters, usually only a few processes contribute, other subprocesses 
are suppressed by the Boltzmann factor. 

The micrOMEGAs package just after compilation does not contain the code for matrix 
elements. They are generated and linked in runtime when needed. To generate the 
matrix elements we use the CalcHEP program [16] in batch mode [38]. The compiled 
matrix elements are stored as shared libraries in the subdirectory 
sources/2-2/ 

The name of the library created corresponds to the names of initial superparticles. Say, 
the library containing XiXi annihilation processes is omglib_ol_ol . so. 

On the first few calls, micrOMEGAs works slowly because it compiles matrix elements. 
After being compiled once, the code for matrix elements is stored on the disk and is 



accessible for all subsequent calls. Each process is generated and compiled only once. 

In case several jobs are submitted simultaneously, a problem occurs when CalcHEP 
receives a new request to generate a matrix element when it has not completed the previous 
one. We delay the operation of the second program. The warning that CalcHEP is busy 
signals the presence of a LOCK file in the directory 
sources/work/ tmp 

If for some reason this file is not removed after the CalcHEP session, the user should 
remove it. 

The executable file generated by micro_make can be moved and executed in other 
directories. However it will always use and update the matrix elements stored in 
micromegas_l . 3 . 0/sources/2-2 

5.3 Linking with other codes and including micrOMEGAs into 
other packages. 

One can easily add other libraries to the micrOMEGAs package similarly to the imple- 
mentation of Isajet described in Section 5.1. One needs to pass the library name to 
the linker via the EXTLIB variable defined in micro_make, by specifying the complete 
path to the library. One can include the micrOMEGAs package into other C, C++, or 
Fortran projects. The function prototypes for C and C++ projects are stored in the 
sources/micromegas .h file. All the routines of our package as well as Suspect and 
FeynHiggsFast routines are stored in 

sources/micromegas . a 
which in turn needs the functions of 

sources/ decay2 . a 

to calculate the widths. The user must pass to the linker the library that supports dy- 
namic loading. The name of this library depends on the Unix platform. One can find this 
name in the micro_make file, it is assigned to the LDDL environment variable. 

To attach micrOMEGAs to a C or C++ project, the user should make sure that the 
library of Fortran functions are also passed to the linker. In the micro_make file this 
library is described by the LDF variable. 

5.4 Running micrOMEGAs 1 . 3: examples. 

The directory micromegas_l . 3 . contains several examples of main programs. The files 
sugomg . c and sugomg_f . f are main programs for the evaluation of the relic density in 
the mSUGRA scenario. 



. /mi cr o_make sugomg . c 

generates the executable sugomg which needs 5 parameters 

./sugomg <m0> <mhf> <a0> <tb> <sgn> 
The sugomg executable also understands three additional input parameters as m t , m h {m h ), 
a s (M z ). The output contains the SUSY and Higgs mass spectrum, the value of the relic 
density, the relative contributions of different processes to 1/Q as well as the constraints 
mentioned in Section 4.4. The list of necessary parameters are written on the screen when 
sugomg is called without specifying parameters. 

. /mi cr o_make sugomg_f . f 
compiles the corresponding Fortran code. In this case the input parameters are requested 
after launching the program: 

> ./sugomg_f 

Enter mO mhf aO tb sgn 
> 

By default these programs call Suspect for solving the RGE equations. One can easily 
change the RGE code by replacing the suspect SUGRA call by the appropriate one in 
sugomg. c or sugomg_f . f. 

The program s.cycle.c performs the calculation over 10 mSUGRA test points [39]. 
Results for these points for all RGE programs mentioned in our paper are presented in 
the file data/s_cycle . res. 

Finally the omg . c and omg_f . f programs evaluated the relic density in the case of 
the unconstrained MSSM. The input parameters are read from a text file written in the 
format of the ewsblnitFile routine. In the C-version the file should be passed as a 
parameter, for example 

./omg data/data03 

If several sets of parameters are passed to the program, the calculation will be done in 
a cycle. The Fortran version also works in a cycle, waiting for a file name as input and 
finishes after an empty line input. 

The directory data contains 22 "data*" test input files for this routine. These param- 
eter sets were chosen to check the program in special difficult cases where either strong 
co-annihilation and/or Higgs pole contribute significantly in relic density evaluation. Re- 
sults of relic density calculation for all these 22 test points using the option when all 
masses are evaluated at tree-level are stored in file data\omg.res. 



5.5 CalcHEP interactive session. 



The CalcHEP [16] program used for matrix element generation is included in the micrOMEGAs 
package. The user can calculate interactively various cross sections both in the general 
MSSM and in SUGRA models. To realize this option the user has to move to the calchep 
subdirectory and launch 
. /calchep 

The implementation of the MSSM and SUGRA models in CalcHEP is identical to the 
one in micrOMEGAs described in previous sections. There are two auxiliary parameters, 
LCOn and dMbOn which switch ON/OFF loop corrections to the MSSM particle spectrum 
and SUSY-QCD correction to h,H,A^ bb decays respectively. If LC0n>0 or dMbOn>0 the 
corresponding correction is taken into account. 

The list of parameters contains also the scale parameter Q which should be set de- 
pending of the scale of the process under consideration. This parameter contributes to 
the running of a s and to the running masses of t and b quarks. Here we use the stan- 
dard MS formulae without including the higher order QCD corrections 5 presented in 
Section 2.1. 

For the SUGRA model, all four RGE packages presented in micrOMEGAs can be used, 
Suspect is defined by default. External RGE packages are available for CalcHEP if they 
were already properly installed in the micrOMEGAs package as described in section 5.1. 
To include another RGE package one has to edit the model in CalcHEP (in the Edit 
model menu). The suspect SUGRA call should be commented in the Constraints menu 
while the line corresponding to the call for another routine should be uncommented. 
The symbol for comment is %. In the Edit model menu one can also defined the non- 
universal SUGRA model. By default, mSUGRA boundary conditions are implemented. 
To modify this, first comment the lines in the Constraints menu which express the GUT 
scale parameters in Table 3 in terms of the mSUGRA parameters. The corresponding 
non-universal parameters should then be introduced as new variables in the Variables 
menu. 

In this realization of MSSM/SUGRA all widths of super-partners are evaluated auto- 
matically at tree-level including all l->2 decay modes generated in the model. The relic 
density and other constrains mentioned in section 4.4 are included in the list of Constrains 
and automatically attached to CalcHEP numerical sessions. 

In CalcHEP numerical sessions for 2->2 processes we provide an option to construct 
a plot for the v ■ a dependence on the incoming momentum. This option is found under 
5 These corrections can be simulated by decreasing of scale Q. 



the Simpson menu function. 
5.6 Sample output file 

Running micrOMEGAsl . 3 with the default values of the standard parameters and choosing 
the Suspect RGE package with the mSUGRA input parameters 

sugomg 107 600 5 1 

will produce the following output: 



Higgs masses and widths 



h 


: Mh 


= 116 





(wh 


=2 


5E-03) 


H 


: MHH 


= 899 


2 


(wHh 


= 1 


9E+00) 


H3 


: MH3 


= 898 


5 


(wH3 


=2 


2E+00) 


H+ 


: MHc 


= 902 





(wHc 


=2 


3E+00) 



Masses of SuperParticles : 



~ol 


MNE1 = 


249 


1 1 


I "11 


MS11 = 


254 


2 I 


I ~eR 


MSeR = 


256 





~mR 


MSmR = 


256 


I 


1 ~nl 


MSnl = 


413 


1 | 


I ~ne 


MSne = 


413 


4 


~nm 


MSnm = 


413 


4 I 


I ~eL 


MSeL = 


420 


2 I 


I ~mL 


MSmL = 


420 


2 


~12 


MS12 = 


420 


4 I 


I "1+ 


MCI 


468 


3 I 


I ~o2 


MNE2 = 


468 


5 


~o3 


MNE3 = 


780 


I 


I "2+ 


MC2 


793 


2 I 


I ~o4 


MNE4 = 


794 


3 


~tl 


MStl = 


946 


7 I 


I ~bl 


MSbl = 


= 1153 


1 | 


I ~b2 


MSb2 = 


= 1187 


8 


~dR 


MSdR = 


1188 


4 I 


1 ~sR 


MSsR = 


= 1188 


4 I 


I ~t2 


MSt2 = 


= 1190 


6 


~uR 


MSuR = 


1194 


8 I 


I ~cR 


MScR = 


= 1194 


8 I 


I ~uL 


MSuL = 


= 1248 


2 


~cL 


MScL = 


1248 


2 I 


I ~dL 


MSdL = 


= 1250 


5 I 


I ~sL 


MSsL = 


= 1250 


5 


~g 


MSG 


1358 


1 1 



















Xf=2.67e+01 0mega=8 . 87e-02 



Channels which contribute to 1/ (omega) more than 1%. 
Relative contrubutions in % are displyed 



a 


~ol 


~ol 


-> 


1 


L 


3% 


~ol 


"11 


-> 


Z 


1 


12% 


~ol 


"11 


-> 


A 


1 


2% 


~ol 


~eR 


-> 


Z 


e 


8% 


~ol 


~eR 


-> 


A 


e 


2% 


~ol 


~mR 


-> 


Z 


m 


8% 


~ol 


~mR 


-> 


A 


m 


11% 


"11 


"11 


-> 


1 


1 


2% 


"11 


"LI 


-> 


A 


Z 


3% 


"11 


"LI 


-> 


A 


A 


8% 


~eR 


"11 


-> 


e 


1 


6% 


~eR 


~eR 


-> 


e 


e 


1% 


~eR 


~ER 


-> 


A 


Z 


2% 


~eR 


~ER 


-> 


A 


A 


6% 


~eR 


~mR 


-> 


e 


m 


8% 


~mR 


"11 


-> 


m 


1 


6% 


~mR 


~mR 


-> 


m 


m 



1% ~mR ~MR -> A Z 
2% ~mR ~MR -> A A 
deltartho=9 . 1 1E-06 
gmuon=3. 12E-10 
bsgnlo=3.85E-04 
bsmumu=3 . 13E-09 
MassLimits OK 



Under the same conditions and for the same set of parameters, running the cross- 
section and branching ratios routines 

cs_br 

will produce the following output: 

Example of some cross sections and widths calculation 
for mSUGRA point m0=107.0,mhf=600.0,a0=0.0,tb=5.0 



Z partial widths 



b 


B 


- 3. 


, 684E- 


01 


GeV 


d 


D 


- 3. 


.703E- 


01 


GeV 


u 


U 


- 2. 


.873E- 


■01 


GeV 


c 


C 


- 2. 


.873E- 


01 


GeV 


s 


s 


- 3. 


,703E- 


■01 


GeV 


1 


L 


- 8. 


.378E- 


■02 


GeV 


nl 


Nl 


- 1 


.670E- 


01 


GeV 


nm 


Nm 


- 1 


.670E- 


•01 


GeV 


ne 


Ne 


- 1 


.670E- 


•01 


GeV 


m 


M 


- 8. 


.397E- 


02 


GeV 


e 


E 


- 8. 


.397E- 


•02 


GeV 



Total 2.436E+00 GeV 

h partial widths 
b B - 2.460E-03 GeV 
1 L - 2.552E-04 GeV 
Total 2.716E-03 GeV 

Cross sections at Pcm=500.0 GeV 
e,E->~l+,~l- 

e,E->~l+(468) ,~l-(468) is 7 . 135E-03 pb 
e,E->~ol,~o2 

e,E->~ol(249) ,~o2(468) is 1 . 130E-02 pb 



6 Results 



We have compared the results obtained with micrOMEGAsl . 3 and those obtained with 
DarkSUSY4.0 for 10 benchmarks mSUGRA points[39]. For this check, we have used 
Isajet7.69, m p t ole = 174.3GeV, a s (M z ) = .1172, m b (m b ) = 4.23GeV. The latter is only 
relevant for the calculation of the Higgs widths. 

As seen in Table 8, the two programs agree at the 3% level except at large tan/3. 
This discrepancy is due to a difference in the width of the pseudoscalar. We recover good 
agreement with DarkSUSY (below 3%) if we substitute their value for the pseudoscalar 
width. 



Table 8: Comparison between micrOMEGAsl . 3 and DarkSUSY4.0 



name 


M 


M 1/2 


A 


tan (3 


sgn(n) 


micrOMEGAsl. 3 


DarkSUSY4.0 


A 


107 


600 





5 


1 


0.0944 


0.0929 


B 


57 


250 





10 


1 


0.124 


0.121 


C 


80 


400 





10 


-1 


0.117 


0.115 


D 


101 


525 





20 


1 


0.0876 


0.0864 


G 


113 


375 





20 


1 


0.133 


0.129 


H 


244 


935 





20 


1 


0.166 


0.163 


I 


181 


350 





35 


1 


0.142 


0.132 


J 


299 


750 





35 


1 


0.102 


0.0975 


K 


1001 


1300 





46 


-1 


0.0893 


0.0870 


L 


303 


450 





47 


1 


0.114 


0.0982 



7 Conclusion 

micrOMEGAsl .3 solves with an accuracy at the percent level the evolution equation for 
the density of supersymmetric particles and calculates the relic density of dark matter. 
All possible channels for annihilation and coannihilations are included and all matrix 
elements are calculated exactly in an improved tree-level approximation that uses pole 
masses and loop-corrected mixing matrices for supersymmetric particles. Loop corrections 
to the masses of Higgs particles and to the partial widths of the Higgs (QCD and SUSY) 
are implemented. These higher-order corrections are essential since the annihilation cross- 
section can be very sensitive to the mass of the particles that contribute to the various 
annihilation processes, in particular near a resonance or in regions of parameter space 
where coannihilations occur. Furthermore, both these processes are often the dominant 
ones in physically interesting supersymmetric models, that is in models where the relic 
density is below the WMAP upper limit. 



The relic density can be calculated starting from a set of MSSM parameters defined at 
the weak scale or at the GUT scale. We provide an interface to the four major codes that 
calculate the supersymmetric spectrum using renormalization group equations. Within 
the context of the mSUGRA model, there are still large uncertainties in the computation 
of the supersymmetric spectrum [42], this of course will have a strong impact on the 
prediction for the relic density [14]. An accurate prediction of the relic density within 
SUGRA models then presupposes a precise knowledge of the supersymmetric spectrum. 

New features of the package also include the computation of cross-sections and decay 
widths for any process in the MSSM with two-body final states as well as an improved 
NLO calculation of the b — > S7 branching ratio and a new routine for the B s — > 
decay rate. 
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Appendix 

A List of functions 

A.l micrOMEGAs functions in C. 

int assignVal (char * name, double val) 
void assignValW(char * name, double val) 
int readLesH(char *fname) 
int ewsblnitFile (char * f name, int LC)} 

int ewsbMSSM (tb , MG1 , MG2 , MG3 , Am , Al , At , Ab , MH3 ,mu , Mil , M12 , M13 , Mr 1 , Mr2 , Mr3 , Mql , Mq2 , Mq3 , 

Mul,Mu2,Mu3,Mdl,Md2,Md3,LC) ; int LC; all other parameters are 'double' 
int xxxxxSUGRA ( tb , MG 1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , Ml 1 , M13 , Mr 1 , Mr 3 , Mql , Mq3 , 
Mul,Mu3,Mdl,Md3) 



All parameters are 'double', 'xxxx' is 'suspect' , 'isajet' , 'softSusy' ,or 'spheno' 



int calcDep(int dMbOn) 

int findVal(char * name, double * val) 

double f indValW(char * name) 

void print Var (FILE *f, int N) 

void printMasses (FILE * f , int sort) 

char * lsp(void) 

double lspmass_() 

double darkOmega (double *Xf,int fast, double Beps) 
double darkOmegaFO (double *Xf,int fast, double Beps) 

double print Channels (double Xf, double cut, double Beps, int prcnt, FILE *f ) 

double deltarho_ (void) 

double bsgnlo_ (void) 

double bsmumu_ (void) 

double gmuon_(void) 

int masslimits_ (void) 

double MbRun (double Q) 

double MtRun (double Q) 

double MbEff (double Q) 

double MtEff (double Q) 

double deltaMb(void) 

double decay2 (char*pln, int k, char*pOutl, char*p0ut2) 
void* newProcess (char* procName, char*libName) 

int inf or22(void*address , int nsub, char*plnl,char*pln2,char*p0utl,char*p0ut2, 

double*ml , double*m2 , double*m3 , double*m4) 
double cs22 (void*address , int nsub, double Pern, double cl, double c2, int*err) 
double annihilation (double v, int k, char * pOutl, char p0ut2) 

micrOMEGAsfunctions in Fortran. 

INTEGER FUNCTION assignVal (name , val) 
SUBROUTINE assignValW(name , val) 

INTEGER FUNCTION readLesH(f name) 
INTEGER FUNCTION ewsblnitFile (f name , LC) 

INTEGER FUNCTION ewsbMSSM (tb , MG1 , MG2 , MG3 , Am , Al , At , Ab , MH3 , mu , Mil , M12 , M13 , 
Mr 1 , Mr2 , Mr3 , Mql , Mq2 , Mq3 , Mul , Mu2 , Mu3 , Mdl , Md2 , Md3 , LC) 

INTEGER FUNCTION xxxxSUGRA (tb , MG1 , MG2 , MG3 , Al , At , Ab , sgn , MHu , MHd , 
Mll,M13,Mrl,Mr3,Mql,Mq3,Mul,Mu3,Mdl,Md3) 

All parameters are 'double'. 'xxxx J is 'suspect' , 'isajet' , 'softSusy' ,or 'spheno' 

INTEGER FUNCTION calcDep (dMbOn) 

INTEGER FUNCTION f indVal (name , val) 

REAL* 8 FUNCTION f indValW(name) 

SUBROUTINE printVar(n) 

SUBROUTINE printMasses (sort) 

SUBROUTINE LSP(name) 

REAL* 8 FUNCTION IspMassO 

REAL* 8 FUNCTION darkOmega(Xf , fast , Beps) 

REAL* 8 FUNCTION darkOmegaFO (Xf , fast , Beps) 



REAL*8 FUNCTION printChannels (Xf , cut , Beps ,prcnt) 

REAL*8 FUNCTION deltarhoO 

REAL*8 FUNCTION bsgnloO 

REAL*8 FUNCTION bsmumuO 

REAL*8 FUNCTION gmuonQ 

INTEGER FUNCTION MassLimits () 

REAL*8 FUNCTION MbRun(Q) 

REAL*8 FUNCTION MtRun(Q) 

REAL*8 FUNCTION MbEff (Q) 

REAL*8 FUNCTION MtEff(Q) 

REAL*8 FUNCTION deltaMbO 

REAL*8 FUNCTION decay2(pln, k, p0utl,p0ut2) 
SUBROUTINE newProcess (procName , libName , address) 

INTEGER FUNCTION infor22 (address, nsub, plnl ,pln2 ,pOutl ,p0ut2 ,ml ,m2 ,m3 ,m4) 
REAL*8 FUNCTION cs22 (address , nsub, Pern, cl, c2 , ERR) 
REAL*8 FUNCTION annihilation(v,k, p0ut,p0ut2) 

The types of the parameters are: 

CHARACTER pln*(*) ,plnl*(*) ,pln2*(*) ,pOutl*(*) ,p0ut2*(*) , 

> name*(*) , fname*(*) , procName ,* (*) , libName (*, *) 
REAL*8 val , Xf , Beps , cut , Pcm , c 1 , c2 , v , Q , ml , m2 , m3 , m4 
REAL* 8 tb,MGl,MG2,MG3,Am,Al,At,Ab,MH3,mu,Mll,M12,M13, 

> Mrl,Mr2,Mr3,Mql,Mq2,Mq3, Mul ,Mu2 ,Mu3 ,Mdl ,Md2 ,Md3 
INTEGER n,k, sort, prcnt, ERR, LC,dmbOn, fast, address [2] 

B Implementation of B(B — > 57) in micrOMEGAs 

The calculation for B(B — > S7) in the MSSM is quite involved and requires that one goes 
beyond one-loop. Most of what is described below, as implemented in micrOMEGAs, is 
in fact just a, unified, compendium of different contributions that have appeared in the 
literature. There is no claim of originality, most expressions are taken verbatim. However 
care has been taken in carefully checking all formulae that have appeared in the literature. 
This has helped, for example, identify a few misprints and typos and allowed to generalise 
some results. By giving the detail of the implementation, it is possible to easily modify 
this routine of the micrOMEGAs code in order to include future new contributions both to 
the SM and the MSSM. Note that we redefine in this routine many parameters used in 
micrOMEGAsl .3, for example the running quark masses, this routine can then be used as 
a stand-alone routine. 



B.l General set-up: From Mw to QED corrections 

Our implementation of the Standard Model contribution follows the work of Kagan and 
Neubert [43] very closely. We however include the effect of a running c quark mass 
heuristically so that our results take into account the latest calculations of Gambino and 



Misiak[44] who advocate the use of the MS charm mass, m c (m b ). The (relevant) operator 
basis is 

7 = ^s L v^b R , 

Os = Y^s L a^t a b R . (B.13) 



which defines 



AC 

H eS = V*V tb ]T CM)OM) ■ (B.14) 



The renormalisation scale fi b in (Eq. B.14) is of order m h and is usually let to vary in 
the range (nib/2, 2rrib). The default value in the code is mj. Varying \i b is one measure of 
the theoretical error. The branching fraction writes 



2 

K NLO (5) x B(B^X c eu) . (B.15) 
By default we take 

B(B -> X c eu) = 0.1045 (B.16) 

The kinematical function, f(zo), is defined as 

f(zo) = l-8z + 8z% -z%- 12z*\nz « 0.542 - 2.23(v^ - °- 29 ) ( B - 17 ) 

with z = (m c /m b ) 2 defined in terms of the pole masses, giving a value in the range, y/z = 
0.29 ± 0.02. For the radiative photon we take a = 1/137.036. 

The factor i^NLo(^) involves the photon energy cut-off parameter 8 that shows up at 
the NLO. In micrOMEGAs this value is set to 5 = 0.9 as is generally assumed in order to 
describe the "total" (fake) branching ratio. With the formulae given below the code can 
be modified in a very straightforward way to take into account the full S dependence. 
-^nlo(^) is decomposed in terms of the Wilson coefficients with leading (LO) and next- 
to-leading (NLO) contributions as 



B(B - X sl ) = 



ba 



nf(zn) 



VtsVtb 



Vrh 



K NLO (S) = E ^(^^ReCf^Cf*^) + S(S) 



a 



i,j=2,7,8 

i<3 



2vr 



Re 



(0)*, 



and 



+ S(5) (2Re[^ em V 6 ) C?>(^ b )} - 4 c L m V 6 ) |C?V b )| 2 ) , (B.18) 



CM = Cf V) + ^ Cf V*>) + Cf m V 6 ) + • • • • (B.19) 

47T a s {fi b ) 



The leading-order coefficients at the low scale \i b of order m b are given by 

COl , 1 / _ 12 _6_\ 



cPM = n* (cf\ mw) + III) + E ^> ^ , (B.20) 

where 77 = a s (m w )/a s (j2 b ), and /i® and <2j are known numerical coefficients [46]. 



626126 ^56281 3 1 

hi = ( , , — , ,-0.6494,-0.0380,-0.0186,-0.0057) 

v 272277 51730 7 14' ' ' ' ; 

14 16 6 12 
a t = (—,—,— ,-—,0.4086,-0.4230,-0.8994,0.1456) 

hf ] = (-0.9135,0.0873,-0.0571,0.0209) 
bi = (0.4086,-0.4230,-0.8994,0.1456) (B.21) 

For the running of a s between the scale M w and /i b we use the SM running with 5 
flavours which, to a very good precision, can be implemented as: 

a s (M z ) ( 116 a s (M z )\n(v s (fi)) \ 23 a s (M z ) 

= V ~ — <iV ) Vsifl) = 1 " HMz/fI) 

(B.22) 

The value of a s (M z ) is read in by the main code micrOMEGAs. For the numerical 
values that we quote in this note, we take the default a s (M z ) = 0.1185. 

The next-to-leading Wilson coefficient at fi b , C^\ii b ) is implemented according to 

37 39 i 



Cl*"b b ) = V^C^(M W ) + I^-^)C^(M W ) 

Z297664 is 7164416 u 256868 37 6698884 sA^™,,, n 

+ 7723 7723 -| 7123 7723 Cv (M W ) 

V 14283 ' 357075 ' 14283 ' 357075 ' / 8 V W) 



+ W ( V * " + TX^E(x) + f l + gw)rf\ (B.23) 



( 



4661194 



8516 



o, 



o, 



-1.9043, -0.1008, 0.1216, 0.0183) 



816831 ' 2217' 

fi = ( -17.3023, 8.5027, 4.5508, 0.7519, 2.0040, 0.7476, -0.5385, 0.0914) 
gi = ( 14.8088, -10.8090, -0.8740, 0.4218, -2.9347, 0.3971, 0.1600, 0.0225) 



and 



„, , :r(18- lLx-rr 2 ) x 2 (15 - 16x + 4x 2 ) , 2, 
E{X) = 12(1 - X )> + 6(1-^ lnX - 3 ]nX - 



(B.24) 



The QED coefficients C\ eTa \fi b ) and k^\fi b ) are 
The result for C? em ^ (fi b ) is 

/-((em)/ \ /32 __9. 40 _1_ 88 16 \ (o) 

C 7 (fib) = [t^V 23 --^V 23 + V 23 ) C 7 J (m w ) 



V75 '' 69 
32 9 
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rj 2 3 + 
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H r? 23 — 

4347 ' 



?7 23 + 



640 



775 



704 



7723 ) C^ 0) (miy) 



1449 ' 1449 ' 1725 

359 _i7 4276 _« 350531 _ji 

77 23 -J- 77 23 -I 77 23 

3105 ' 121095 ' 1009125 ' 
5956 _6_ 38380 m 748 ie 

77 23 -4- fj 23 — 77 23 



15525 



169533 



8625 



(B.25) 



4rV 6 ) = ^ fa- 1 - 1) 



23 V ' ' 7T /i b 

The coefficient functions kij(8, /if,) in (B.18) are given by [43] 



(B.26) 



k 77 (5, fi b ) 



^ J77(C>) +6(0) — , 



7T 



2rr 



^27(5, /ife) 
fc 78 (5, // 6 ) 



5(5) 
5(5) 



/i 6 / 9 m^ 



«s(/fb) 

2tt v 



/(«>) 

a s (/4>) 



- 1 



6A 2 



+ 



7T 



fl b J VT 



Oisiflb) 



7T 



{i,j} = {2,2}, {8,8}, {2,8}, 



with the Sudakov factor 



S{5) = exp 



2 -^(ln 2 5+ 7 -ln5 

3vr V 2 



(B.27) 



(B.28) 



and 777 = y, 727 = and 737 = — ^ are entries of the anomalous dimension matrix. 
The value of the hadronic parameter is A 2 = 0.12GeV 2 . 

For these functions we deviate slightly from KN[43] in the sense that we define Zq 
in terms of the pole masses in the kinematics factor and also z that differs from z by 
the use of the MS running charm mass, m c (mj), as advocated recently by Gambino and 
Misiak[44] in order to reduce the NNLO uncertainty. We then take 

v^ = 0.22±.04 (B.29) 

everywhere else in Eq. B.27. 

The other coefficients are given by 



r 7 



10 8tt 2 

T 



Re(r 8 ) 



44 8tt 2 



9 ' v 7 9 27 ' 

Re(r 2 ) w -4.987+ 12. 78(v^- 0.22), R(z) « 3.672 - A.U(y/z - 0.22)(B.30) 

Note that the scale jib, in Eq- B.27, of relevance in semileptonic B decays is in principle 
different from the one in the radiative decay. 

The real-gluon radiation functions fij(5) can be coded for a general photon energy 
cut-off. They are taken from [43] and write as 



frr(S) 
fss(5) 



25 3 

105 + 5 2 - — + 5(5-4) In 5 
2tt 2 




1 

3 L 

^ |4Li 2 (l - <J) - ^- + 8m(l-<5) -5(2 + 5) In 5 

25 3 1 
+ 75 + 35 2 - — - 2 [25 + 5 2 + 41n(l - 5)] ln bs \ , 

In — we take — ~ 50 , 

m s m s 

Li 2 (l -5) 5\n5-\ 1 

v ; 6 4 4 12 



x \zj 2 



Jdx (1 — x)(l — #,$) 


3/ 28 (5) = -y fdx(l-x 5 )Re 



; t<4, 



2 In 



-2arctan 2 v /t/(4-t) 
(v^+v^)/^-^) 2 ;f>4. 



;b.si) 



;b.32) 



Since we will specialise to the case 5 = 0.9, it is more efficient to quote the correspond- 
ing values of , and give approximations to f 2 7, /22 that we use in the code: 



/ 77 (0.9) = 3.20599, 

/ 88 (0.9) = 1.31742, 

/ 78 (0.9) = 0.387341, 

/ 22 (0.9) « 0.107636 - 0.208484^ - 0.156146,2 = 0.05421 - 0.2772e z - 0.156146c 2 , 

/ 27 (0.9) = -3/ 28 (0.9) « -0.190805 + 0.948865^ - 0.787805,2 



-0.02023 + .6020e z - 0.7878e: 



e z = v^-0.22 



(B.33) 



Taking z = 0.22 2 , Eq.B.29, rather than z = 0.29 2 mainly affects K 27 especially through 
Re(r 2 ). Note the large coefficient of e z in Re(r 2 ) in Eq. B.30. 



B.2 Standard Model contribution 

This contribution we take from [47]. We first define the functions 



x(7-5x-8x 2 ) x 2 (3x-2), 

+ — tt" -TT- lux 



24(x - l) 3 

x{2 + 5x — x 1 
8(x- l) 3 



4(x - l) 4 

3x 2 
A(x- l) 4 



lnx 



(B.34) 
(B.35) 



B.2.1 LO at M w 

We have 

Cff M W) = F 7 ] 8 (x tw ) (B.36) 
where x tlu is defined in terms of the running top mass at the weak scale, /iw — M w . 

m 2 t (n w ) 



M 2 W 



;b.37) 



For the NLO top-quark running mass at the scale jiw we follow [47] 



m t {^w) = rn t (m t ) 
m t (m t ) = m t 



a s (fi w ) 


23 


_a s (m t ) _ 




Aa s {m t ) 




3 71 


J 



ajmj 7o m /TT _ Pi\ ( <*s(nw) _ ^ 
4tt 2/3 1,70™ #J ^.(m,) ',. 
a s (m t ) 



m t (m t ) 2 = m 2 



1 



3 TT 



(B.38) 



m t is the pole mass which in this note we take as m t = 174.3 ± 5.1GeV for comparison 
with other authors. 

Po - y , Pi - — , To - 8 > 7i - — ^— (B.39) 

B.2.2 NLO SM 

We have 

Cg SM (M w ) = G 7)8 (x tw ) (B.40) 



-16x 4 - 122x 3 + 80x 2 - 8x T . /„ 1 \ 6x 4 + 46a: 3 - 28x 2 , , 



G ' W = L H 1 "J + 3(x-l)» WX 

-102x 5 - 588x 4 - 2262x 3 + 3244x 2 - 1364a; + 208 , 

ma; 



81(x — l) 5 

1646a; 4 + 12205a; 3 - 10740a; 2 + 2509a; - 436 ,„ . 

+ m^W (R41) 



„, s -4x 4 + 40x 3 + 41x 2 + x T . / 1\ -17a; 3 - 31a; 2 , , 

G 8 (x) = -. r-; Li 2 1 M -. -z — In 2 x 

8V ; 6(x- l) 4 V xJ 2(x- l) 5 

-210a; 5 + 1086a; 4 + 4893a; 3 + 2857a; 2 - 1994a; + 280 , 

H ; lux 

216(x- l) 5 

737a; 4 - 14102x 3 - 28209a; 2 + 610a; - 508 ,„ , 

+ 7 B.42 

1296(x-l) 4 V ; 

B.2.3 Results and Comparisons 

To check the different components of the SM part, we have also introduced the Bij func- 
tions [43]. These can be very useful if one wants to introduce the effects of New Physics 
through the Wilson coefficients defined at the scale M w . Dismissing any right-handed 
light quark operator and assuming purely real contributions, the New Physics contribution 
can be written as 

CiT = x 7>s C™> SM (M W ) (B.43) 

then the are the coefficients of the different Xj factors (linear and quadratic, and a Xj 
independent term). In other words, the contribution of the New Physics can be expressed 
as 



= .622 + B 27 x 7 + B 28 x 8 + .677 x 2 + B 88 Xg + B 78 x 7 x 8 

(BAA) 



Note the assumption in [43] that the proportionality factor is the same for the LO and 
NLO. In our case we define a larger set of By by allowing 

C?f = x% C?^ SM (B.45) 

As a check on the SM results note that we exactly recover the results for the Wilson 
coefficients at Mw for both 6*7,8 and at both LO and NLO (this also agrees with [47]. More 
satisfying is that we recover all the results of [43]. Below, see Table 9, we show our results 
for the coefficients B^. Here we switch to the default values of [43] (\/z = m c /m b = 0.29, 
m b = 4.80 GeV, m t = 175 GeV , a s (m z ) = 0.118, \V*V tb \/\V cb \ = 0.976) and fi b = fx b 

Table 9: Values of the coefficients Bij(6) in units of 10 -4 , for different choices of fi b 



li b 5 


B 2 2 B 77 B 88 B 2 7 B 2 s B 78 




m b /2 0.90 
0.30 
0.15 


1.322 0.335 0.015 1.265 0.179 0.074 
1.169 0.322 0.005 1.196 0.136 0.070 
1.081 0.309 0.004 1.144 0.126 0.067 


3.190 
2.898 
2.730 


m b 0.90 
0.30 
0.15 


1.258 0.382 0.015 1.395 0.161 0.083 
1.239 0.361 0.005 1.387 0.137 0.080 
1.200 0.347 0.004 1.354 0.132 0.077 


3.293 
3.210 
3.114 


2m b 0.90 
0.30 
0.15 


1.023 0.428 0.015 1.517 0.132 0.092 
1.041 0.402 0.004 1.552 0.118 0.091 
1.021 0.386 0.004 1.534 0.115 0.088 


3.206 
3.209 
3.149 



Switching back to our default central values, with m b = 4.8GeV, |^*Vtf,|/|Vd,| 2 = 0.971, 
m c = 1.25GeV, m t = 174.3, a fl (Mf ) = 0.1185 and with fi b = Ji b = m b = ml 3 = 4.80GeV 
we find 

10 4 b^ lo (5 = 0.9, yfz = 0.22) = 1.512 + 1.417 a;? + 0.155 x° 8 

+ 0.136 x\ + 0.017 4 + 0.283 (x° 7 f+ 0.014 (x° 8 f + 0.064 x° 7 x° 8 
+ 0.103^x 7 + 0M3x° 7 xl + 0.007^4 + 0.001^4 (B.46) 

whereas in the assumption of [43] we get 

10 4 B* LO (5 = 0.9, yfz = 0.22) = 1.512 + 1.553 x 7 + 0.173 x 8 + 0.386 x 2 7 

+ 0.015 4 + 0.084 x 7 x 8 (B.47) 

leading to 

B NLO,SM( S = ^ = 0>22 ) = 3 723 1Q -4 ( BAg ^ 



This result agrees at the 2 per-mil with the more sophisticated analysis of [44] and is in 
very good agreement with the current experimental value. 



B.3 Charged Higgs contribution 

The charged Higgs contribution in our code is estimated at the M w scale and is based 
on [47], therefore we neglect any small running from M H ± to M w . Indeed either M H ± is 
very large in which case the Higgs contribution is too small, or there is little running in 
a region where a s is not too large. 



B.3.1 LO 

We have 

with 
and 



B.3.2 NLO 

We have 



Ml 



C\ 1)H± {^w) = G"{x Ht ) + A?(x Ht ) In - ^E H (x Ht ) 



Ml 9' 



H 

2 



M 2 H 6" 



(B.49) 



x m = , (B.50) 



to\, . x(3 — 5x) x(3x — 2) , 

^ = 4^-27^^ (B.51) 



Ci 1)H± (»w) = Gf(x m ) + A»(x Ht )\n^-lE H (x Ht ) (B.52) 



G?{x) = 



4 

--x 
3 



4(-3 + 7x-2x 2 ) T . / IN 8-14x-3x 2 l2 



2(-3-a;+12a; 2 -2a; 3 ) 



lna; + 



xj 3(x-l) 4 
7- 13a: + 2x 2 ' 



1 2 



3(x-l) 4 

x(18 - 37x + 8x 2 



(x - iy 



■Lio 



tan/5 2 9' [x 

-50 + 251a; - 174a; 2 - 192a; 3 + 21a; 4 , 

H -. r-= mx 

9{x-l) 5 



1\ x(-14 + 23a; + 3a; 2 ) , 2 

1 - - + — ? —z - In 2 a; 

x ) \x — \y 



+ 



797 - 5436a; + 7569a; 2 - 1202a; 3 



108(rr- l) 4 



(B.53) 



2 

--x 
9 



21 - 47a; + 8x 2 2(-8 + 14a; + 3a; 2 ) , 

- H ; — In x 



(x - iy 



(x - iy 



1 2 



tan/3 2 9 



x 



-31 - 18a; + 135a; 2 - 14a; 3 x(U - 23a; - 3a; 2 ) , 

1 ; — In x 



6(x-l) 



(x-iy 



B.54) 



Gf(x) 



1 



-36 + 25a; - 17a; 2 ) / 1\ 

2(x- iy 12 v ~ x ) + (x - iy 



19 + 17a;, 2 
In x 



-3 - 187a; + 12a; 2 - 14a; 3 ) , 3(143 - 44a; + 29a; 2 ) 
H 77 77 In x + — 



4(x - l) 4 



x-r 



l l 



tan (3 2 6 



a; 



x{30- 17a; + 13x 2 ) / I 
— L12 ( 1 



(x - iy 



X 



x(31 + 17a;) , 2 

~7 TT^ ln x 

[x — 1 r 



-226 + 817a; + 1353x 2 + 318a; 3 + 42x 4 , 

H ; — ma; 

36(a;- l) 5 

1130 - 18153a; + 7650a; 2 - 4451a; 3 

+ 216(x- l) 4 



(B.55) 



Af(x) = --x 



81 - 16a; + 7a; 2 19 + 17a; 



2(x- iy 



{x-iy 



\nx 



1 1 
-x 



E H {x) 



tan/5 2 6 
1 



-38 - 261a; + 18a; 2 -7a; 3 x(31 + 17x) 



6fx- l) 4 



+ 



(a; - 1) E 



In a; 



tan/5 2 



a;(16 - 29a; + 7a; 2 ) x(3x - 2) , 

+ 777 7T7- 



36(a; - l) 3 



6(x- 1) 



(B.56) 
(B.57) 



As a check we recover exactly the LO and NLO values quoted in Table 1 of [47]. 

B.4 SUSY contributions 

We only consider the contribution from charginos (and accompanying squarks). Here we 
follow [45] rather closely but adapt the expressions for the e&,e;,(i). 



B.4.1 LO 



We first consider the LO SUSY contribution at the SUSY scale. Though the code allows 
to choose this scale, our default value is set at the gluino mass /U susy = nig. This is also 
the scale we take for the SUSY Am b corrections. 



With 



Fi 3 \x) 



F<?\x) 



5 — 7x x(3x — 2) , 

+ 777 7TT In X, 



6(x-l) 2 3(a;-l) 3 

1 + X X 



2(x-iy {x-iy 



In a; 



(B.58) 
(B.59) 



J2M, 2 

2 / - 

77 C t 'Kl - StK2 f- . . 2 

3 \ v2sm/3MvF / 



W 7^(1) 



tl 



3 \ V 2 sin /3 M w / m~ 



+ 



1 / U a2 V al M w 



cos (3 \ y/2: 



m 



.2 F (3) 



2 irO), 



7,8 x 



Ua2 Va2 m t 

2 sin P m 



(3), 



(B.60) 



with obvious notations for the sparticles. 1/m 2 here and in the following. Our 

diagonalising matrices for the chargino is as in [45] ( as well as our convention for the sign 
of n) 

M 2 



M w y/2smp\ 

\M w y/2cos/3 



The squark mixing and definitions are also as in [45] 
Cq = cos 9g , Sq = sin 9g with 

qi = Cqq L + Sq q R , q 2 = -Sq q L + Cqq R , and > rriq 2 



(B.61) 



(B.62) 



B.4.2 Am b corrections and large tan/? effects 

m b = V2M W — cos (3 (1 + e b tan/3) , 5m h = 

We implement in our code as follows. First define 
a; In a; y In y 



Anib 



e^tan /?. 



where (x,y)ij = m 2 /m 2 . 



(B.63) 



H(i,j,k) = H 2 (x ik ,y jk ) (B.64) 



g(m 1 ,m 2i g 2 ) = 7j + + ( /^ )2 - M^'/Q 2 )) , ^m'M (B.65) 



Then 



uKVsusy) ST TT V ~ At tal1 ^ Hff + T> 
+ ^a 2 2^ U a2 M (t 1 ,t 2 ,X a ) V a2 

16 71 a=i,2 rn xt 
+ a ^ Z \ M 2 tan/5 \-%-H(M 2 , fi, h) + \h{M 2 , fi, t 2 ) 

2 2 \ 

+ ^rH(M 2 , fi, h) + ^rH(M 2 , h, b 2 ) (B.66) 

bi b 2 ) 

Note that we have included the electroweak contributions with enhanced tan/5. For 
this part we make the approximation that the masses of the charginos are given by M 2 
and fi, and neglect the mixing matrices, we have also neglected the U(l) contribution 
which is even smaller. Although the formulae above include non-leading tan/3 effects, 
for consistency we have not coded these contributions in our program. The electroweak 
corrections agree with those in [49] whereas the sbottom term is missing in [50] . However 
we find the SU(2) gauge contribution to be rather small compared to the strong and 
Yukawa contributions. 



We turn now to e' b (t) 

2a s (fi susy ) A b j tan /3 - fi 
3 7r m 



= + ~ c i c b H{t l ,b 2 ,g) + c i s- b H{t l ,b l ,g) + 



9 

s 2 l clH{t 2 ,b 2 ,~g) + s 2 slH{t 2 ,b l ,. 



. Uti^susy) v-^ , r * A t /i/tan/3 r 2 2 uf r 7 0\ , 2 2 tt (7 7 0\ , 
V ~^~A L N oA C t C b H{t 2 , h, Xa) + C~ t S~ b H(t 2 , b 2 , Xa) + 



10^ ^ " m x o 

2 2 Tj/7 t n\ . 2 2 



stclHit^xl) + stslH^b^xl) 



N, 



a .3 



2 



2 2 \ 

+ ^2 H ^ *o + 4k H{M2 ' /i? * 2) (R67) 

Note that we have added an electroweak contribution, in the same approximation as 
in 5m b . Most importantly, the sign of the Yukawa contribution and the elements of the 
diagonalising matrices that appear in Eq. B.67 are different from those in [45]. We have 
verified this by explicit calculation of e' b (t), moreover with our formula, we find that in 
the decoupling limit we do indeed have e' b (t) — > e b , which would not have been the case 
had we blindly used the formula of [45]. One of the authors of [45], Paolo Gambino, has 
recently confirmed our implementation. 



As for e' t (s), we find that only the QCD contribution remains. Note, as said previously, 
we shall only keep the tan (3 enhanced term in our code. 



For e' t (s) 



2a s /j, + A t / tan/3 



3tt 



c~H(t 2 ,s,g) + s~H(t u s,g) 



(B.68) 



We find no tan (3 enhanced electroweak gauge contribution. 

Once more let us stress that, although we have derived, for the e, e' the tan (3 enhanced 
and the non tan f3 enhanced we will only keep the tan f3 enhanced terms, since these are 
the ones that can be resummed. Therefore in our numerical analysis that takes this 
resummation into account we only keep the tan (3 enhanced terms. 

As we mentioned earlier these e's contributions are to be evaluated at the scale Q 2 = 
fi 2 usy > m 2 which we associate with the gluino mass. In particular a s is to be evaluated 
here taking into account 6 active quarks. 

-l 

(B.69) 



Vs = Oi s (nsusY)/a s (iJ,w) = ( 1 + l^LJ^L \n(jj, SUSY /fiw) 



2tt 

Also the Yukawa (and top mass) that is used for the chargino contribution is 

Vtinsusv) = yt(fJ>w) 
x , 



ois(vsusy) 


4/7 


a s (m t ) 


12/23 


a s (m t ) 




a 8 (n w )_ 





'1 + 



9j/i (rat) 
87ra s (m t ) 



as(HSUSY) 
a s (rn t ) 



1/7 



(B.70) 



At all scales we relate y t to m t as follows 

yKQ 2 ) = 



, . 2 . _ 2na{M z ) 1 + tan (3 2 m 2 (Q 2 ) 



s 2 



tan (3 



2 M 2 
m w 



(B.71) 



We used a fixed tan/3 at all scales, we neglect the small running of a between Mz and 

f^susy 

Though the e effects are to be extracted at /i susy , they are included in the SM and 
charged Higgs Wilson coefficient at /iw 

8Cj S ^ (leading tan f3)(/j, w ) ■ 



h-e' b (t)}tanf3 (2) 

r 7,8 \ x tw) 



SC^ ^ (leading tan (3){nw) 



l + tan j3 
[e' t (s) + e b ] tan/3 (2) 

7— 7 a ^7,8 K x Ht)- 

1 + €b tan (3 



[B.72) 
^B.73) 



For the chargino contribution we first evaluate the Wilson coefficients at the scale 
A*susy by including the top mass effect and the e effects. In our implementation we 



assume, as occurs in most cases in mSUGRA, that both stops are heavy, so that they are 
decoupled together at the same scale /xsusy, 



+ 



a=l,2 



3 m 2 



2 / ~ ~ m t (vsusY 

- \CfVai - S~ t V a2 



2 Ml 



V2sin(3M w ) m? *W x *xt> 



, 2 M, 2 



- i:\StVai + ciV a2 —7=——— . 2 
3 \ V 2 sin /3 M w / m| 



ti 

2 

W 7^(1) 



K b U a2 V al M w 



cos/3 I y/2: 



rn 



.2 F (3) 



2 17(3)/ 



U a 2 V a2 mtinsusr) 



t w 



2 sin /? m 



(3), 



Xa 



*2 Xa- 



^B.74) 



with 



X 6 = 1/(1 + e 6 tan /3) 



(B.75) 



These are then evolved to /xjy as 



C7 (mw) 



16 

' 3/3' /nrX 



(B.76) 



At 

Hw ah the contributions (SM, H ± including tan/5 effects) are added together with 
those of the chargino contribution and evolved according to Eq. B.20. 

To check on the SUSY part and the implementation of tan (3 enhanced terms, we have 
first verified that we had perfect agreement with Fig. 4 of [45]. This in fact is only a check 
on the implementation of the e's in the SM and H + contribution with fixed e. A full check 
of the SUSY contribution requires a quite large set of inputs. In [45] one can read the 
effect on B — > X s ^ of a SUGRA model. Since the outputs of SUGRA codes can differ 
quite a bit, we have not used our own RGE but requested the weak scale parameters from 
Paolo Gambino used for their Fig. 2 and Fig. 3 in [45] , which in passing includes only 
the a s contribution to the e. We have found an excellent agreement both for positive and 
negative //. 
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